Human CTB BPS Project

Overall Objective

Load Libraries

library(tidyverse)
library(cowplot)
library(broom)
library(plotly)

Import data

Convert data from ‘wide’ to ‘long’ format

# data
data1 <- data %>%
  gather(Sample,Count,2:145)

# Separate samples by identifiers 
data2 <- data1 %>% 
  separate(Sample, into=c("Sample_ID","Dilution_factor",
                          "Injection","Tech_rep", sep = "_")) %>% 
  select(-`_`)

# Standards
standards1 <- standards %>% 
  gather(Sample,Count,2:49)


standards2 <- standards1 %>% 
  separate(Sample, into=c("Sample_ID","When","Dilution_factor",
                          "Nano_day","Injection","Tech_Rep", sep = "_")) %>% 
  select(-`_`)

Factor the data into categorical variables

# Refactoring Columns for samples
data2$Sample_ID <- as.factor(data2$Sample_ID)
data2$Dilution_factor <- as.numeric(data2$Dilution_factor)
data2$Injection<- as.factor(data2$Injection)
data2$Tech_rep <- as.numeric(data2$Tech_rep)


# Refactoring COlumns for key
key$Sample_ID <- as.factor(key$Sample_ID)
key$Time <- as.factor(key$Time)
key$Treatment <- as.factor(key$Treatment)
key$Volume <- as.numeric(key$Volume)
key$Patient_ID <- as.factor(key$Patient_ID)

key$Treatment <- factor(key$Treatment,levels = c('DMSO','EGF','BPS','BPS_EGF'))
key$Patient_ID <- factor(key$Patient_ID,levels = c('1','5','7'))


# Refactoring columns for standards
standards2$Sample_ID <- as.factor(standards2$Sample_ID)
standards2$When <- as.factor(standards2$When)
standards2$Dilution_factor <- as.numeric(standards2$Dilution_factor)
standards2$Injection <- as.factor(standards2$Injection)
standards2$Nano_day <- as.numeric(standards2$Nano_day)

Back calculate standards

standards2 <- standards2 %>% 
  mutate(True_Count=Dilution_factor*Count)

# Set the correct order of 'categorical factors'
standards2$Nano_day <-  factor(standards2$Nano_day, levels=c('1','2','3','4'))
standards2$When <- factor(standards2$When, levels=c('before','after'))

Summarize three technical standard replicates

standards3 <- standards2 %>% 
  group_by(particle_size,Sample_ID,When,Dilution_factor,Nano_day,Injection) %>% 
  summarise( tech_N = length(True_Count),
             tech_mean = mean(True_Count),
             tech_sd = sd(True_Count),
             tech_se = tech_sd/sqrt(tech_N))
standards3

Summarize standards by injection

standards4 <- standards3 %>% 
  group_by(Nano_day,When,particle_size) %>% 
  summarise( inj_N = length(tech_mean),
             inj_mean = mean(tech_mean),
             inj_sd = sd(tech_mean),
             inj_se = inj_sd/sqrt(inj_N))
standards4

Plot before and after plots, facet by experimental day

std_plot <- standards4 %>% 
  ggplot(aes(x = particle_size, y = inj_mean, color=When))+
  geom_line(size=2) + xlim(0,300)+ #line size, x-axis scale
  geom_ribbon(aes(ymin=inj_mean-inj_se, ymax=inj_mean+inj_se),
              alpha=0.4,fill = alpha('grey12', 0.2)) + #error bars
  scale_y_continuous(expand=c(0,0))+ #set bottom of graph
  xlab("Particle Size") + # X axis label
  ylab("\nMean Particle Concentration\n(particles/ml)\n") + # Y axis label
  ggtitle("Nanosight Histogram of\n100nm Standards\n(Facets Represent Experimental Days)")+ #title
  labs(color="Condition")+ #Label table title
  facet_wrap(~ Nano_day)

std_plot

## Warning: Removed 1400 rows containing missing values (geom_path).

# ggsave("Standards_histogram_plot.png",
#        height = 5, width = 7, dpi = 300, units= "in")

Standards particle concentrations from each experimental day

standards_df <- standards4 %>% 
  group_by(Nano_day,When) %>% 
  summarise(total=sum(inj_mean))

standards_df

Bar graph of standards particle concentrations

standards_bar <- standards_df %>% 
  ggplot(aes(x=Nano_day,y=total,fill=When))+
  geom_col(position="dodge")+
  scale_y_continuous(expand=c(0,0))+ #set bottom of graph
  xlab("Experimental Day") + # X axis label
  ylab("\nMean Particle Concentration\n(particles/ml)\n") + # Y axis label
  ggtitle("Nanosight Histogram of\n100nm Standards")+ #title
  labs(color="When") #Label table title

standards_bar

# ggsave("Standards_bar_plot.png",
#        height = 5, width = 7, dpi = 300, units= "in")

Intraassay variability

Intra.assay_cv <- standards_df %>% 
  group_by(Nano_day) %>% 
  summarise(Intra_Day_N = length(total),
            Intra_Day_mean = mean(total),
            Intra_Day_sd = sd(total),
            Intra_Day_se = Intra_Day_sd/sqrt(Intra_Day_N),
            Intra_Day_cv = Intra_Day_sd/Intra_Day_mean )
Intra.assay_cv

# # Save as .csv
# write_csv(Intra.assay_cv,"Intra.assay_cv.csv")

Interassay variability

Inter.assay_cv <- Intra.assay_cv %>% 
  summarise(Inter_Day_N = length(Intra_Day_mean),
            Inter_Day_mean = mean(Intra_Day_mean),
            Inter_Day_sd = sd(Intra_Day_mean),
            Inter_Day_se = Inter_Day_sd/sqrt(Inter_Day_N),
            Inter_Day_cv = Inter_Day_sd/Inter_Day_mean )
Inter.assay_cv

# # Save as .csv
# write_csv(Inter.assay_cv,"Inter.assay_cv.csv")

Sample analysis

Back calculate the original concentration of the sample

data2 <- data2 %>% 
  mutate(True_Count = Dilution_factor*Count)
data2

Average three technical readings

data3 <- data2 %>% 
  group_by(particle_size,Sample_ID,Dilution_factor,Injection) %>% 
  summarise( tech_N = length(True_Count),
             tech_mean = mean(True_Count),
             tech_sd = sd(True_Count),
             tech_se = tech_sd/sqrt(tech_N))
data3

Summarize samples by injection (average both injections)

data4 <- data3 %>% 
  group_by(particle_size,Sample_ID,Dilution_factor) %>% 
  summarise( inj_N = length(tech_mean),
             inj_mean = mean(tech_mean),
             inj_sd = sd(tech_mean),
             inj_se = inj_sd/sqrt(inj_N))
data4

# Average technical replicates and merge with key
merge <- left_join(key,data3, by= "Sample_ID")

merge

# Save as .csv
# write_csv(merge,"Technical_replicate_average.csv")
 
# Average injection replicates and merge with key
merge1 <- left_join(key,data4, by= "Sample_ID")

merge1

# #Save as .csv
# write_csv(merge1,"Injection_replicate_average.csv")

Quick visualizations

Graphing all samples at 48 hours

sample_plot_48 <- merge %>%
  filter(Time == "48") %>% 
  ggplot(aes(x=particle_size, y=tech_mean,color=Injection ))+ #plot
  geom_ribbon(aes(ymin=tech_mean-tech_se,
                  ymax=tech_mean+tech_se),
                  alpha=0.2,fill = alpha('grey12', 0.2)) + #error bars
  geom_line(size=2.0, alpha = 0.8) + xlim(0,500)+ #line size, x-axis scale
  scale_y_continuous(expand=c(0,0))+ #set bottom of graph
  xlab("Particle Size") + # X axis label
  ylab("\nMean Particle Concentration/ml\n") + # Y axis label
  ggtitle("Nanosight Histogram of\nhCTBS treated with BPS (48 hours)")+ #title
  labs(color="Injection")+ #Label table title
  facet_grid(Patient_ID ~ Treatment)+
  geom_vline(xintercept = 200)+
  annotate("text", x= 350, y = 1E8, label= "200nm")

sample_plot_48

# ggsave("Nanosight_Sample_Histogram_48hr.png", plot = sample_plot_48,
#        height = 10, width = 14, dpi = 200, units= "in")

Graphing all samples at 96 hours

sample_plot_96 <- merge %>%
  filter(Time == "96") %>% 
  ggplot(aes(x=particle_size, y=tech_mean,color=Injection ))+ #plot
  geom_ribbon(aes(ymin=tech_mean-tech_se,
                  ymax=tech_mean+tech_se),
                  alpha=0.2,fill = alpha('grey12', 0.2)) + #error bars
  geom_line(size=2.0, alpha = 0.8) + xlim(0,500)+ #line size, x-axis scale
  scale_y_continuous(expand=c(0,0))+ #set bottom of graph
  xlab("Particle Size") + # X axis label
  ylab("\nMean Particle Concentration/ml\n") + # Y axis label
  ggtitle("Nanosight Histogram of\nhCTBS treated with BPS (96 hours)")+ #title
  labs(color="Injection")+ #Label table title
  facet_grid(Patient_ID ~ Treatment)+
  geom_vline(xintercept = 200)+
  annotate("text", x= 350, y = 5E7, label= "200nm")

sample_plot_96

## Warning: Removed 1000 rows containing missing values (geom_path).

# ggsave("Nanosight_Sample_Histogram_96hr.png", plot = sample_plot_96,
#        height = 10, width = 14, dpi = 200, units= "in")

Particle concentration values for each of the samples

merge2 <- merge1 %>% 
  group_by(Time, Treatment, Volume,Patient_ID) %>% 
  summarise(particle_conc=sum(inj_mean))
merge2

Correct for resuspension volume

merge3 <- merge2 %>% 
  mutate(particle_count = (Volume/1000)*particle_conc, # Create new column with number of particles
         corrected_particle_conc = (particle_conc/1E9)) # Create new column with correct particle concentration
merge3

# Save as .csv
# write_csv(merge3,"Adjusted_particle_concentration.csv")

Barplot

merge3$Patient_ID <- factor(merge3$Patient_ID,levels = c('1','5','7'))
merge3$Time <- factor(merge3$Time,levels= c('48','96'))

plot1 <- merge3 %>%
  group_by(Time,Patient_ID,Treatment) %>% 
  ggplot(aes(x = Treatment, y = corrected_particle_conc, fill = Time, group=Time)) +
  geom_bar(aes(text = paste("Particle Concentration:",
                            corrected_particle_conc)),
            position = "dodge", stat= "identity")+
  xlab("\nTreatment\n") + # X axis label
  ylab("\nMean Vessicle Concentration\n(10^9 particles/ ml)\n") + # Y axis label
  ggtitle("Effect of BPS on Extracellular Vessicle\nRelease of hCTBs\n")+
  scale_y_continuous(breaks = seq(0,14,2),
                     limits = c(0,14),
                     expand = c(0,0))+ # set bottom of graph
  labs(fill="Time(hr)")+ # Label table title
  facet_wrap(~Patient_ID)

## Warning: Ignoring unknown aesthetics: text

plot1

# ggsave("BPS_treated_hCTBs_sample_facet_plot.png",
#       height = 8, width = 10, dpi = 300, units= "in")

Double Facet Plot

sample_double_facet <- merge3 %>%
  group_by(Time,Patient_ID,Treatment) %>% 
  ggplot(aes(x = Treatment, y = corrected_particle_conc, fill = Time, group=Time)) +
  geom_bar(aes(text = paste("Particle Concentration:",
                            corrected_particle_conc)),
            position = "dodge", stat= "identity")+
  xlab("\nTreatment\n") + # X axis label
  ylab("\nMean Vessicle Concentration\n(10^9 particles/ ml)\n") + # Y axis label
  ggtitle("Effect of BPS on Extracellular Vessicle\nRelease of hCTBs \n")+
  scale_y_continuous(breaks = seq(0,14,2),
                     limits = c(0,14),
                     expand = c(0,0))+ # set bottom of graph
  labs(fill="Time(hr)")+ # Label table title
  facet_wrap(Time~Patient_ID)

## Warning: Ignoring unknown aesthetics: text

sample_double_facet

# ggsave("BPS_treated_hCTBs_sample_double_facet_plot.png",
#       height = 8, width = 10, dpi = 300, units= "in")

Boxplot

plot2 <- merge3 %>%
  group_by(Time, Treatment) %>% 
  ggplot(aes(x = Treatment, y = corrected_particle_conc, color = Patient_ID )) +
  geom_boxplot(color = "black", fill = NA)+
  geom_point(aes(text = paste ("Patient_ID",Patient_ID)),position = 'jitter', size = 5)+
  xlab("\nTreatment\n") + # X axis label
  ylab("\nMean Vessicle Concentration\n(10^9 particles/ ml)\n") + # Y axis label
  ggtitle("Effect of BPS on Extracellular Vessicle\nRelease of hCTBs\n")+
  labs(fill= "Time (hr)")+
  facet_wrap(~Time)

## Warning: Ignoring unknown aesthetics: text

plot2

 # ggsave("BPS_treated_hCTBs_48_96_plot.png",
 #       height = 10, width = 14, dpi = 500, units= "in")

Interactive Plot

ggplotly(plot2)

merge4 <- merge3 %>% 
  group_by(Time,Treatment) %>% 
  summarise(   N = length(corrected_particle_conc),
            mean = mean(corrected_particle_conc),
              sd = sd(corrected_particle_conc),
              se = sd/sqrt(N))

merge4

averaged_plot <- merge4 %>% 
  ggplot(aes(x = Treatment, y = mean, fill = Time )) +
  geom_bar( position = "dodge", stat= "identity")+
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=.5, 
                size=0.8, colour="black", position=position_dodge(.9)) + #error bars
  scale_y_continuous(breaks = seq(0,12,3),
                     limits = c(0,12),
                     expand=c(0,0))+ #set bottom of graph
  xlab("\nTreatment\n") + # X axis label
  ylab("\nMean Vessicle Concentration\n(10^9 particles/ ml)\n") + # Y axis label
  ggtitle("Effect of BPS on Extracellular Vessicle\nRelease of hCTBs (n = 3)\n")+
  labs(fill= "Time (hr)")+
  facet_wrap(~Time)

averaged_plot

# ggsave("Averaged_BPS_treated_hCTBs_48_96_plot.png",
#       height = 5, width = 7, dpi = 500, units= "in")

Statistics

tidy(shapiro.test(merge3$particle_conc))

fit <- aov(corrected_particle_conc ~ Time * Treatment, data=merge3)

tidy(fit)

tukey_fit <- TukeyHSD(fit)

tukey <- tidy(tukey_fit)
tukey

tukey %>% 
  filter(adj.p.value <0.05) %>% 
  arrange(adj.p.value)

Filtering on 200nm or less

nano_data <- data3 %>% 
  filter(particle_size <200.5) %>% 
  group_by(particle_size,Sample_ID,Dilution_factor) %>% 
  summarise( inj_N = length(tech_mean),
             inj_mean = mean(tech_mean),
             inj_sd = sd(tech_mean),
             inj_se = inj_sd/sqrt(inj_N))
nano_data

# Average technical replicates and merge with key
merge_nano <- left_join(key,nano_data, by= "Sample_ID")

## Warning in left_join_impl(x, y, by$x, by$y, suffix$x, suffix$y): joining
## factors with different levels, coercing to character vector

merge_nano

merge_nano2 <- merge_nano %>%  
  group_by(Time, Treatment, Volume,Patient_ID) %>% 
  summarise(particle_conc=sum(inj_mean))

merge_nano3 <- merge_nano2 %>% 
  mutate(particle_count = (Volume/1000)*particle_conc, # Create new column with number of particles
         corrected_particle_conc = (particle_conc/1E9)) # Create new column with correct particle concentration
merge_nano3

nano_plot <- merge_nano3 %>% 
 group_by(Time,Patient_ID,Treatment) %>% 
  ggplot(aes(x = Treatment, y = corrected_particle_conc, fill = Time, group=Time)) +
  geom_bar(aes(text = paste("Particle Concentration:",
                            corrected_particle_conc)),
            position = "dodge", stat= "identity")+
  xlab("\nTreatment\n") + # X axis label
  ylab("\nMean Vessicle Concentration\n(10^9 particles/ ml)\n") + # Y axis label
  ggtitle("Effect of BPS on Extracellular Vessicle (<200nm)\nRelease of hCTBs \n")+
  scale_y_continuous(breaks = seq(0,14,2),
                     limits = c(0,14),
                     expand = c(0,0))+ # set bottom of graph
  labs(color="Condition")+ # Label table title
  facet_wrap(~Patient_ID)

## Warning: Ignoring unknown aesthetics: text

nano_plot